Graded mesh B-spline collocation method for two parameters singularly perturbed boundary value problems

The solutions of two parameters singularly perturbed boundary value problems typically exhibit two boundary layers. Because of the presence of these layers standard numerical methods fail to give accurate approximations. This paper introduces a numerical treatment of a class of two parameters singularly perturbed boundary value problems whose solution exhibits boundary layer phenomena. A graded mesh is considered to resolve the boundary layers and collocation method with cubic B-splines on the graded mesh is proposed and analyzed. The proposed method leads to a tri-diagonal linear system of equations. The stability and parameters uniform convergence of the present method are examined. To verify the theoretical estimates and efficiency of the method several known test problems in the literature are considered. Comparisons to some existing results are made to show the better efficiency of the proposed method. Summing up:• The present method is found to be stable and parameters uniform convergent and the numerical results support the theoretical findings.• Experimental results show that the present method approximates the solution very well and has a rate of convergence of order two in the maximum norm.• Experimental results show that cubic B-spline collocation method on graded mesh is more efficient than cubic B-spline collocation method on Shishkin mesh and some other existing methods in the literature.


Introduction
Differential equations having small perturbation parameters multiplying the highest and higher derivative terms are called singularly perturbed problems (SPPs).The solutions of these types of problems exhibit multi scale characters.That is, there are a narrow region called boundary layer in which their solution changes rapidly and the outer region where solution changes smoothly.Classical numerical methods are known to be inadequate to solve such multi scale problems unless extremely fine meshes are used.However, this fine uniform mesh requirement by most of the classical methods for stability in the singularly perturbed case leads significantly large algebraic system.Hence this incorporates the massive computational cost.Therefore, design of robust layer adapted meshes for such problems is an important task.In this paper the following problem is considered: with boundary conditions where  0 ,  1 are given scalars and the two small perturbation parameters 0 < ,  ≪ 1 are such that  2  → 0 as  → 0 (  2 ≪ ) or   2 → 0 as  → 0 (  ≪  2 ).The functions  (  ) ,  (  ) , and  (  ) are assumed to be sufficiently smooth functions and satisfying Such types of problems are called two parameters singularly perturbed boundary value problems (TPSPBVPs).The applications of such types of problems appear in transport phenomena arising in chemistry and biology [1] , in lubrication theory [2] , in chemical reactor theory [3] , and so on.
When  = 0 , the problem (1) -( 2) reduces to the well-studied reaction diffusion problem of single parameter singularly perturbed boundary value problem (SPSPBVP) while for  = 1 it becomes the convection diffusion problem.For the problem classes of these two special cases, there is a huge number of research papers devoted to their numerical solving, however there are only few number of research papers for the purely two parameters case 0 < ,  ≪ 1 .For problem classes of case  = 0 (reaction diffusion problem)for instance research papers like [4][5][6][7][8] and for  = 1 (convection diffusion problem) papers like [9][10][11][12] are available in the literature.
The solution of problem ( 1) -( 2) exhibits narrow regions called inner(boundary layer) regions where it behaves irregularly and changes rapidly and the outer region where it behaves regularly and changes slowly.Because of the presence of the boundary layers the classical numerical methods are known to be unsuitable & fail to give accurate results when the perturbation parameters are small unless extremely fine meshes are used.But this fine mesh requirement by the classical methods leads to a very large algebraic system which incorporates high computational cost.Therefore it is important to develop suitable parameters uniform layer adapted numerical methods to get approximate solutions of such problems.Parameters uniform numerical methods are methods which satisfy error bounds of the form where  , is exact solution of the continuous problem,   is numerical approximation of the exact solution  , , || ⋅ || is the maximum pointwise norm, N is the number of mesh intervals used and c is a positive constant independent of N and the parameters.A numerical method is said to be parameters uniformly convergent of order p if For further details, the reader is referred to the books and research papers available in the literature.Among the books [13][14][15][16][17] are mentioned.Some of the research papers available in the already existing literature are expressed as follows.O'Malley in [18] studied the nature of the two parameter problem asymptotically and examined that the solution behavior depends on the ratios   2 as  → 0 and  2  as  → 0 .Roos & Uzelac in [19] constructed a streamline diffusion finite element method on Shishkin mesh and shown that the method is almost second order parameters uniformly convergent.The authors in [20] proposed upwind difference method on Shishkin-Bakhvalov mesh to solve a TPSPBVP and found the method to be first order parameters uniform convergent.Linß & Roos in [21] used a simple upwind difference scheme on Shishkin mesh to solve a TPSPBVP and found the method to be parameters uniformly convergent of almost order one.The authors in [22] have proposed B-spline collocation method on Shishkin mesh & shown that it has a second order parameters uniform convergence.The Ritz-Galerkin method has been proposed by Kadalbajoo et.al in [23] .These authors have used a piecewise uniform Shishkin mesh to resolve the boundary layers and the method has been found to be almost second order parameters uniform convergent.
The authors in [24] presented a comparative study of fitted mesh finite difference, Ritz-Galerkin finite element and B-spline collocation methods on Shishkin mesh.They have shown that B-spline collocation method has second order parameters uniform convergence, the Ritz-Galerkin finite element method has almost second order parameters uniform convergence and the finite difference method has almost first order parameters uniform convergence.In [25] a semi-linear TPSPBVPs is solved using exponential spline on a Shishkin mesh and the method is shown to be parameters uniformly convergent.The authors in [26] developed a quadratic Bspline collocation method for TPSPBVPs on exponentially graded mesh and they have shown the method to be parameters uniformly convergent of order two.
Besides the fact that there are fewer number of research papers for TPSPBVPs as compared to the number of research papers on SPSPBVPs, most of the research papers on TPSPBVPs are on Shishkin meshes.As far as our knowledge of the literature the number of research papers available on graded meshes is very few as compared to those on Shishkin mesh.Moreover, for some methods like finite difference method on Shishkin mesh it is known that rate of convergence is deteriorated by the presence of logarithmic factor whereas graded mesh methods which are without logarithmic factors give better rate of convergence.Hence the aim of this paper is to develop and analyze a graded mesh cubic B-spline collocation method using a spline interpolant which satisfies the same boundary conditions for approximating solutions of TPSPBVPs.Experimental results in this paper shown that cubic B-spline collocation method on graded mesh (present method) is better than cubic B-spline collocation method on Shishkin method (existing method).The flow chart for the work is given in the appendix.

Continuous problem solution properties
In this section bounds on the solution and its derivatives of problem (1) -(2) which will be used in error analysis will be established.Characteristic equation of the homogeneous equation corresponding to the differential equation in (1) is The solutions of Eq. ( 6) become These two real solutions  0 (  ) < 0 &  1 (  ) > 0 describe the boundary layers at x = 0 & x = 1 respectively.Define Recall that the solution behavior of problem ( 1) -( 2) depends on the ratios  2  → 0 as  → 0 and   2 → 0 as  → 0 .
case (1):  2  → 0 as  → 0 .In this case one can easily verify that  0 =  1 = √   , which gives: Hence behavior of solution in this case is similar to that of the reaction diffusion problem (  = 0 ) and exhibits boundary layers at x = 0 & x = 1 with boundary widths of (  − 1 2 ) [27,28] for each layer.case (2):   2 → 0 as  → 0 .For this case it is possible to show that Since the parameters  and  are non-zero small numbers and using Eq.(10) it can be observed that −1 =  2   0 implies: Hence when   2 → 0 as  → 0 , the layer at x = 1 is stronger than the layer at x = 0 and the boundary layer widths are (  −1 ) and (   ) at x = 0 and x = 1 respectively [27,28] .Remark 1.For case (2) when  (  ) < 0 , a strong layer at x = 0 exists when compared to the layer at  = 1 .(Itmeans when  (  ) < 0 ,  1 ≪  0 ).Moreover for both case (1) and case (2) it can clearly be shown that The operator L in Eq. ( 1) satisfies the following maximum principle.
The following stability estimate is obtained.

Mesh construction
This section presents the graded mesh as constructed in [26] .Define the constants  0 and  1 by where  ∈ (0 , 1) &  ≥ 1 is a parameter to be chosen by the user.Let  = {   ∶   =   ,  = 0 , 1 , ⋯ , } be a uniform mesh.Let us construct a non-uniform mesh Ω = {   }   =0 , where N is the number of sub-intervals in [0,1] which is a multiple of 4,  ∈ ℕ &  ≥ 8 .Let the two mesh generating functions Φ 0 & Φ 1 which are important to resolve the boundary layers be defined by Then the mesh points   ,  = 0 , 1 , ⋯ ,  are defined as: Subdividing the interval [0,1] into three subintervals: . Let the mesh step size ℎ  be defined by ℎ  =   −   −1 ,  = 1 , 2 , ⋯ , .The mesh step size ℎ  for  = 1 , 2 , ⋯ ,  satisfies the following properties (which will be used in error analysis): Defining   =  −Φ  ,  = 0 , 1 ; it can be proved that the mesh widths in the layer regions satisfy: Cubic B-splines on non-uniform mesh

B-splines of degree 0
The B-splines of degree 0 denoted by  0  are defined by B-splines of degree k as in [30] are defined by the recurrence relation Then the recurrence relation (24) above can be re-written in the form: Since  0  is a piecewise polynomial of degree zero and since    is linear,  1  is a piecewise polynomial of degree less than or equal to one.Similarly through same reasoning,     is a piecewise polynomial of degree less than or equal to k.
This means the support for This implies that B-splines are positive and never negative.
For cubic B-splines k = 3, the recursive formula ( 24) above becomes: and the support for  3  (  ) is (   ,   +4 ) .By recursively using (27) for  3  , it is possible to obtain: Following B-spline values which will be used latter in the next section can be obtained from the above definition of B-splines.
The numerical scheme In this section, a cubic B-spline collocation method on the graded mesh is derived.

𝐴𝛾 = 𝑑
where  is an (N+1) by (N+1) tridiagonal nonsingular matrix,  = (  0 ,  1 , ⋯ ,   )  is the vector of the unknowns &  is the vector of the right hand side terms.The elements of the tridiagonal matrix  = (   ) are obtained from the coefficients of  0 and  1 of Eq. ( 46) , from the coefficients of  −1 ,   and  +1 of Eq. (47) for  = 1 , 2 , ⋯ ,  − 1 and from the coefficients of  −1 and   of Eq. ( 48) .The elements of  that are found this way satisfy: Hence the coefficient matrix  is strictly diagonally dominant and therefore nonsingular. is nonsingular in turn implies that the system  =  can be solved for  0 ,  1 , ⋯ ,   and use Eqs.( 44) and (45) to calculate for  −1 and  +1 .Thus the cubic B-spline collocation method has a unique solution on problem (1) -( 2) which is given by Eq. (35) .

Convergence analysis
In this section error bound for cubic B-spline collocation method is found.Let the mesh intervals be denoted by   = [   −1 ,   ] ,  = 1 , 2 , ⋯ , .For ,  ∈ ℕ ,  <  , let: , where Π  is the space of polynomials of degree less than or equal to p.
the sum of all  ′   equals 1 by Lemma (5) .
Using theorem (50) and above inequality and this depicts that the method is theoretically parameters uniform convergent of order two. □

Numerical examples
In this section the proposed numerical scheme is applied on a few test problems with two small parameters having two boundary layers to validate the theoretical results. .
The exact solution for this equation is , where . For Examples (1) and ( 2) for any value of N the maximum pointwise absolute errors   , and the order of convergence   , are calculated by: where  , and   , are the exact solution and the numerical solution of the problems in the examples.The exact solution for Example (3) is not available.As the exact solution is not available for Example (3) , the accuracy of its numerical solution will be computed using the double mesh principle.This principle is defined as follows: For any fixed value of , the maximum point-wise error   , of the numerical solution  (   ) will be calculated by where    is the numerical solution at the node   of the original graded mesh with  number of mesh intervals and  2  2  is the numerical solution at the node  2  on a mesh containing 2  number of mesh intervals (i.e. a mesh containing the mesh points of the original graded mesh and its mid points given by

Results & discussion
In the previous section, the proposed cubic B-spline collocation method on graded mesh is applied on three test problems known in the literature to demonstrate its efficiency in obtaining the numerical solution of a class of TPSPBVPs of the form (1) - (2) .The maximum point-wise errors and the orders of convergence for examples (1), ( 2) and (3) for fixed values of , different values of  and for  = {64 , 128 , 256 , 512 , 1024 , 2048} are shown by Tables 1 -6 .It is observed from these results that as  increases the maximum absolute errors decrease independently of the values of the perturbation parameters, which shows the parameters uniform convergence of the method.The log-log plot of the maximum pointwise absolute errors vs N for Example (1) is given in Fig. 8 .In this plot the lines corresponding to the various values of  are overlapping and it seems that the plot has only a single line.As one moves from left to right (as N increases) along this plot the lines are getting down (pointwise maximum absolute errors are decreasing) independently of the values of the perturbation parameters.Further, from the plot it can be seen that the rate of convergence is approximately two (this can also be seen from Tables 1 -6 .Hence the log-log plot also illustrates the parameters uniform convergence of the present method and that its rate of convergence is two.Hence Tables 1 -6 and Fig. 8 validate the theoretical finding that the present method       the error in the present method which is 8 .7025  −04 is smaller than the error 9 .0080  −03 in method [22] , and the rate of convergence 2.0037 in the present method is higher than the rate of convergence 1.6351 in [22] , and so on.Smaller error shows that the method is more accurate and a higher rate of convergence shows that the method is quicker in converging to the exact solution.Hence Table 7 shows that the proposed method is more accurate and having a higher rate of convergence than the method in [22] , that is, the cubic B-spline collocation method on graded mesh (present method) is better (more efficient) than the cubic B-spline collocation method on Shishkin mesh (method in [22] ).Similarly, Table 8 shows that the present method is more accurate and having a higher rate of convergence (more efficient) than the method in [31] already existing in the literature.Although it is possible to observe in a    similar way that the method in [26] is better in accuracy than the present method from Table 9 , it can also be seen from the same table that the present method has higher rate of convergence than the present method.That is, the present method has improved both accuracy and rate of convergence of the methods in [22] and in [31] whereas it improved only the rate of convergence (not accuracy) of the method in [26] .

Conclusion
From the discussion in the above section it can be concluded that, the method is shown to: approximate the solution very well, be parameters uniform convergent of order two in the maximum norm, be more efficient than some methods in the existing literature and adaptable for computer implementation.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.